Numerical simulation is widely used in industrial fields as a method of simulating a physical system by using a computer. In most cases, there is a desire to model the transport processes occurring in the physical system. What is being transported is typically mass, energy, momentum, or some combination thereof. By using numerical simulation, it is possible to model and observe a physical phenomenon and to determine design parameters, without actual laboratory experiments and field tests.
The principle of numerical simulation is to numerically solve equations describing a physical phenomenon by a computer, such as fluid flow. Such equations are generally ordinary differential equations and partial differential equations. These equations are typically solved by linearizing the equations and using numerical methods such as the finite element method, the finite difference method, the finite volume method, and the like. In each of these methods, the physical system to be modeled is divided into smaller gridcells or blocks (a set of which is called a grid or mesh), and the state variables continuously changing in each gridcell are represented by sets of values for each gridcell. In the finite difference method, an original differential equation is replaced by a set of algebraic equations to express the fundamental principles of conservation of mass, energy, and/or momentum within each gridcell and transfer of mass, energy, and/or momentum transfer between gridcells. These equations can number in the millions. Such replacement of continuously changing values by a finite number of values for each gridcell is called “discretization.” In order to analyze a phenomenon changing in time, it is necessary to calculate physical quantities at discrete intervals of time called timesteps, irrespective of the continuously changing conditions as a function of time. Time-dependent modeling of the transport processes proceeds in a sequence of timesteps. Many numerical schemes have been used for solving the partial differential equation, including for example the well known Godunov scheme. In the Godunov scheme, conservation variables are considered as piecewise constant over the mesh cells at each time step and the time evolution is determined by the exact solution of what is known as the Riemann problem. To reduce computation time, the changes in the model parameter values defined for individual grid cells are typically piecewise-constant within the cells rather than continuously varying.
The Wikipedia article “Godunov's Scheme” as accessed on Apr. 9, 2009 (http://en.wikipedia.org/wiki/Godunov%27s_scheme) describes Godunov's scheme as follows:                “In numerical analysis and computation fluid dynamics, Godunov's scheme is a conservation numerical scheme, suggested by S. K. Godunov in 1959, for solving partial differential equations. In this method, the conservative variables are considered as piecewise constant over the mesh cells at each time step and the time evolution is determined by the exact solution of the Riemann problem (shock tube) at the inter-cell boundaries (Hirsch, C. (1990), Numerical Computation of Internal and External Flows, vol 2, Wiley.).        Following Hirsch, the scheme involves three distinct steps to obtain the solution at t=(n+1) Δt from the known solution at t=n Δt, as follows:        Step 1 Define piecewise constant approximation of the solution at t=(n+1) Δt. Since the piecewise constant approximation is an average of the solution over the cell of size Δx, the spatial error is of order Δx, and hence the resulting scheme will be first-order accurate in space. Note that this approximation corresponds to a finite volume method representation whereby the discrete values represent averages of the state variables over the cells. Exact relations for the averaged cell values can be obtained from the integral conservation laws.        Step 2 Obtain the solution for the local Riemann problem at the cell interfaces. This is the only physical step of the whole procedure. The discontinuities at the interfaces are resolved in a superposition of wave satisfying locally the conservation equations. The original Godunov method is based upon the exact solution of the Riemann problems. However, approximate solutions can be applied as an alternative.        Step 3 Average the state variables after a time interval Δt. The state variables obtained after Step 2 are averaged over each cell defining a new piecewise constant approximation resulting from the wave propagation during the time interval Δt. To be consistent, the time interval Δt should be limited such that the waves emanating from an interface do not interact with waves created at the adjacent interfaces. Otherwise the situation inside a cell would be influenced by interacting Riemann problems. This leads to the CFL condition |amax|Δt<Δx/2 where |amax| is the maximum wave speed obtained from the cell eigenvalue(s) of the local Jacobian matrix.        The first and third steps are solely of a numerical nature and can be considered as a projection stage, independent of the second, physical step, the evolution stage. Therefore, they can be modified without influencing the physical input, for instance by replacing the piecewise constant approximation by a piecewise linear variation inside each cell, leading to the definition of second-order space-accurate schemes, such as the MUSCL scheme.”        
In geologic modeling, a digital representation of the detail internal geometry and rock properties of a subsurface earth volume, such as a petroleum reservoir or a sediment-filled basin is built. In the oil and gas industry, geologic models provide geologic input to reservoir performance simulations which are used to select locations for new wells, estimate hydrocarbon reserves, plan reservoir-development strategies, and/or perform other hydrocarbon extraction activities. The spatial distribution of permeability is a key parameter in characterizing reservoir performance, and, together with other rock and fluid properties, determines the producibility of the reservoir. For sandstone reservoirs, the spatial distribution of permeability is a function of the grain-size distribution of sands which compose the reservoir, the compartmentalization of those sands by barriers of finer grained material, and the mineralogy and burial history of the reservoir.
Most geologic models built for petroleum applications are in the form of a three-dimensional array of model blocks (cells), to which geologic and/or geophysical properties such as lithology, porosity, acoustic impedance, permeability, and water saturation are assigned (such properties will be referred to collectively herein as “rock properties”). The entire set of model blocks represents the subsurface earth volume of interest. The goal of the geologic-modeling process is to assign rock properties to each model block in the geologic model.
The geologic modeling process can use many different data types, including but not limited to rock-property data obtained from cores, well logs, seismic data, well test and production data, and structural and stratigraphic surfaces that define distinct zones within the model. Other data types may include surface topographical data and aerial and satellite images that may be used to infer or predict subsurface geologic formations by observing and/or predicting similar geologic activities (erosion, deposition, fluid flow, sediment flow, etc.) on the earth's surface. In general, the resolution or spatial coverage of the available data is not adequate to uniquely determine the rock properties in every geologic model cell. Other assumptions about the variability in these properties are made in order to populate all model cells with reasonable property values. Geocellular techniques, object-based modeling, and process modeling are three main ways to populate the discretized geologic volume with properties.
In the geocellular approach, the relationship between properties in nearby cells is specified statistically. Geostatistical estimation methods (which may be either deterministic or probabilistic) are used to compute rock property values within cells. These methods take into account distance, direction, and spatial continuity of the rock property being modeled. Deterministic estimation methods commonly calculate a minimum-variance estimate of the rock property at each block. Probabilistic estimation methods develop distributions of the rock-property values and produce a suite of geologic model realizations for the rock property being modeled, with each realization theoretically being equally probable. The spatial continuity of a rock property may be captured by a variogram, a well-known technique for quantifying the variability of a rock property as a function of separation distance and direction. U.S. Pat. Nos. 5,838,634, 6,381,543 and 6,480,790 discuss geocellular modeling methods embodied in processing flows which include repetitive optimization steps to drive the geologic model toward conformance with geologic and geophysical data types such as wells, seismic surveys and subsurface fluid production and pressure data. Most commercial modeling software packages, including PETREL, GOCAD and STRATAMODEL, contain a wide spectrum of geostatistical tools designed to fulfill the requirements of reservoir geologists and engineers. While these methods can readily accommodate data control points such as wells and geophysical constraints such as seismic data, they generally do not closely replicate the geologic structures observed in natural systems.
In the object-based approach, subsurface reservoir volumes are treated as assemblages of geologic objects with pre-defined forms such as channels and depositional lobes. U.S. Pat. No. 6,044,328 discloses one object-based modeling scheme that allows geologists and reservoir engineers to select geologic objects from an analog library to best match the reservoir being modeled. The appropriateness of the analog is judged by the operator of the process based on their geologic experience. Most commercial software packages including PETREL, IRAP-RMS and GOCAD implement objects as volumetric elements that mimic channels and lobes using simplified elements based on user deformable shapes such as half pipes and ellipses. Other examples of object-based models are the model depositional objects, such as a river belt and the braided stream network, are placed sequentially on top of each other according to some algorithms. While these models try to mimic the real depositional structures, they do not attempt to capture the physics of water flow and sediment transport that, over geologic time, determined the rock properties at a particular subsurface location.
Process-based models attempt to reproduce subsurface stratigraphy by building sedimentary deposits in chronological order relying to varying degrees on a model or approximation of the physical processes shaping the geology. Although process models are rarely used in current industrial practice, U.S. Pat. Nos. 5,844,799, 6,205,402 and 6,246,963 describe three such methods. These methods employ diffusion or rule-based process models to create basin-scale models with limited spatial detail inadequate for reservoir performance simulation. Another type of model is hydrodynamics-based gridding, which in one aspect may generate grid surfaces conforming to time surfaces during the process of sedimentation. Hydrodynamics-based gridding is the subject of U.S. patent application Ser. No. 11/629,822, “Method for Geologic Modeling Through Hydrodynamics-Based Gridding (Hydro-Grids)”, having inventors Li et al. and filed Dec. 15, 2006.
Other types of process-based models that could potentially provide the most accurate representation of geologic structures are physics-based, process-based models. These models are based on basic physics of fluid flow and sediment transport, and build the subsurface stratigraphy by explicitly modeling the fluid flow, sediment transport, erosion and deposition that occurred during the formation of the reservoir. Compared with other types of geologic models, these physics-based process-based models can better preserve the spatial continuities of sedimentary body geometrics as well as flow baffles and barriers such as shale and mud deposits.
A component in the physics based process-based model is to solve the fluid flow and sediment transport equations. Often this involves solving shallow water or depth-averaged turbidity current equations. A problem in solving shallow-water and depth-averaged turbidity current equations using Godunov-type schemes is the imbalance between flux gradients and bed slope source terms due to the discretized representations of flow and the underlying topography. For stationary flow problems, an ideal variable reconstruction and flux integration scheme should produce a sum of normal components of flux terms evaluated at the edges of a computational cell that exactly balances bed slope source terms evaluated at the center of the cell, and should produce zero changes in flow momentums in the cell. In practice however, these two terms generally do not balance each other in most of the low resolution schemes, as well as in many high resolution schemes. The error introduced by the imbalance is often cumulative, and can significantly degrade the accuracies of the schemes, especially in locations where there are significant topographic variations and controls. An improved method of predicting fluid flow, thereby to improve production management of hydrocarbon resources, is desirable.